Neutral Mutations and Punctuated Equilibrium in Evolving Genetic Networks 
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Boolean networks may be viewed as idealizations of biological genetic networks, where each node 
is represented by an on-off switch which is a function of the binary output from some other nodes. 
We evolve connectivity in a single Boolean network, and demonstrate how the sole requirement of 
sequential matching of attractors may open for an evolution that exhibits punctuated equilibrium. 
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Evolution of life is presumably a random process with 
selection |J. It has been discussed whether this pro- 
cess can be viewed as some hill climbing process [||, 
or whether evolution mostly happens as a random walk 
where changes do not influence the phenotype, and thus 
are neutral || . Originally, the case of evolution as adap- 
tation in an externally imposed fitness landscape has 
been proposed by Sewall Wright ||, and later formed 
the basis for models of punctuated equilibrium by New- 
man (§ and Lande §. The case for neutral evolution has 
been presented by Kimura || , and is experimentally sup- 
ported on the micro-level by the observation that there 
are many functionally identical variants of most of the 
important macromolecules of life. 

The observation of punctuated equilibrium in the fos- 
sil record, recently discussed by Gould and Eldredge 
||, may be taken as an indication that evolution of a 
species consists of exaptations of jumping from one hill 
top to another nearby in some fitness landscape. Nat- 
urally such jumps will be rare, separated by large time 
intervals where species are located at a fitness peak, and 
the resulting evolutionary pattern will show punctuations 
as indeed seen in the fossil record. This picture of sin- 
gle species evolution in a given fixed landscape has been 
modeled explicitly by Newman j| and Lande ||. 

However, also neutral evolution may show punctua- 
tions as, for example, might be visualized by finding the 
exit in a labyrinth or from finding a golf hole by means 
of a random walk in a flat landscape. The picture here 
is that genetic changes always take place, but that the 
phenotypic changes only rarely occur. This has recently 
been demonstrated by modeling the evolution of RNA 
secondary structure by P. Schuster and collaborators . 
For these molecules, mutating a single nucleotide often 
does not induce any changes in their secondary structure, 
and the mutation is considered neutral. Occasionally, 
however, one mutation can lead to a complete readjust- 
ment of the structure, usually accompanied by a major 
change in its functionality. 

In any case, as demonstrated by Bak and Sneppen ||, 
punctuated equilibrium on the organism level might be 
connected to the episodic punctuations observed on the 



ecosystem level. The crucial element of such an extrapo- 
lation is that the environment of each species depends on 
species which are ecological neighbors, thereby allowing 
punctuations to propagate across the ecosystem. 

In the present paper we propose to evolve a single ge- 
netic network, ideally representing a single species. The 
evolution is driven by a noisy environment. The evo- 
lutionary step consists of random mutations combined 
with selection of mutants preserving the phenotype with 
respect to a given environment. Thus, the only require- 
ment in this minimalistic model is continuity in pheno- 
type. Other changes in genotype are allowed, creating a 
path of neutral mutations. We will discuss how this re- 
quirement of continuity in evolution may constrain and 
guide the evolution of an individual species in the face of 
a constantly changing environment. 

Our fundamental constituents are the genes of the or- 
ganism, and the evolution we consider is on the genetic 
network level. Although genetic networks consist of bio- 
chemical switches ||, it has been proposed that the on- 
off nature of these switches can be well approximated by 
Boolean functions jp!o - 12 1 , eventually with asynchronous 
updating Jl3[ |. We here consider networks of random 
Boolean functions, idealized by synchronized updating. 
The functionality we test for is attractors of these net- 
works Boolean networks are known to exhibit a 
rich dynamical behavior, including fixed points, periodic 
attractors and long transients. Further, the number of 
attractors, their length, and the length of the transients 
strongly depend on the connectivity number ]l4|] . In this 
paper we do not address any question about the timescale 
of these attractors. Instead we consider a longer evolu- 
tionary timescale connected to the change in geometry of 
the networks under mutation. 

We implement continuity in evolution by testing for 
reaching a given attractor on subsequent steps, but al- 
lowing changes that modify attractors that are not tested 
from the actual initial condition. In subsequent steps, 
the initial condition (modeling the environment) assumes 
new random values which subsequently allow previous 
neutral mutations to surface in the phenotype. The phi- 
losophy of this subdivision between initial state and func- 
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tion is that while the Boolean functions are manifested 
by various DNA binding proteins, the initial state of the 
system is set by the chemical composition of the envi- 
ronment. This environment changes due to conditions 
beyond the control of the Boolean gene regulatory cir- 
cuit. 

Consider a genetic network with N genes. Each of 
these genes can be assigned a Boolean variable oi = or 
1. For each of the N genes we define an updating ma- 
trix in the form of a lookup table which determines its 
output for each of the possible 2^ input states from the 
N genes in the system. This Boolean updating matrix is 
assigned random values, all rules are a priori equally pos- 
sible. The matrix is effectively quenched on evolutionary 
timescales. Finally we define which gene is actively con- 
nected to which, by a matrix Wij that defines the input 
to gene number i from gene number j as Wij<jj. The en- 
try value of the connectivity matrix can take values 
0, if z is not connected to j, and 1, if i is connected to 
j. Typically only a fraction of the connectivity matrix 
entries is in use, and the average number of connected 
inputs per gene is called the connectivity K. It varies 
between and N, meaning that K may include self cou- 
plings. Thus K — means that all is fixed to the output 
state specified by input (0, ...,0) to all genes. 

The system we evolve is the set of couplings w^j in a 
single Boolean network. The simulation starts with a low 
but finite connectivity, here an initial average connectiv- 
ity of K = 1 per site. One evolutionary time step of the 
network is: 

1. Select a random input state to the network {ci}. It- 
erate the system, called the mother, from this state 
until a final attractor is determined. 

Create a daughter network by a) adding, b) remov- 
ing, or c) adding and removing a weight in the cou- 
pling matrix Wij at random, each option occurring 
with probability p — 1/3. Iterate the daughter sys- 
tem from the same initial state as that selected for 
the mother and test whether it reaches the same 
attractor as the mother system did. In case it does 
then replace mother with daughter network and go 
to step 3. In case another attractor is reached, keep 
mother network and go to step 3. 

Then finally one random bit of the total N * 2 N 
lookup table entries is flipped to another value. 
This allows for a convenient self averaging of the 
system, and in fact represents a very slow change. 

Iterating these steps makes an evolutionary algorithm 
that represents the sole requirement of continuity in 
evolution and how this may proceed under an environ- 
ment that fluctuates. No selective pressure is applied. 
Step 3 rarely affects the active part of the network be- 
cause K typically remains low compared to N, and thus 
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2 K << 2 N . It is included in order to obtain a self av- 
eraging of the system, that elsewise tends to have small 
effective statistics even for the small N we can simulate. 

In Fig. la we show how the connectivity (K) of this 
system evolves with time in a network of size N = 16. 
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FIG. 1. Evolution of the Boolean network connectivity 
with time (a) and closeup on a part of the connectivity evolu- 
tion (b). Note that periods of approximate stasis and sudden 
punctuations appear on both timescales. A single network 
of size N = 16 has been simulated starting with an initial 
average connectivity of K = 1 active inputs per node. The 
connectivity matrix as well as the Boolean updating matrix 
were chosen completely randomly, with all possible Boolean 
rules allowed. Apart from the slow adjustment of Boolean 
rules under step 3 in the model, the system thus evolves in 
a quenched "landscape" of Boolean rules. The connectivity 
K shown is directly measured from the connectivity matrix 
of the network. The effective connectivity [15], defined as 
K minus the number of connected inputs that do not con- 
tribute due to specific Boolean updating matrix entries, has 
somewhat lower values but shows similar overall behavior. 

One observes that the typical K of the network is con- 
fined to lower values than of random networks. This is 
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further quantified in Fig. 2 where the distributions of 
average connectivities are displayed in the statistically 
stationary state. 



show this distribution averaged over the simulation. 




FIG. 2. Distributions of connectivity K in the statistically 
stationary state, obtained from the same simulation as in Fig- 
ure 1. The frequency of connectivities of the new "species" 
is shown ("new K"), as well as the time averaged distribu- 
tion of connectivities of all mutated "species". Note that for 
higher values of if, a mutant network is less likely to match 
the phenotype of its parent. 

Notice that there are two distributions, one counting 
the frequency of connectivities for all new "species" and 
one counting the time averaged distribution. These two 
distributions diverge strongly for high K, because the 
few species with high K have very long lifetimes, i.e., 
for high K it is difficult to find mutations which do not 
change the activity pattern of the networks. In our case, 
the activity pattern consists of the transient and the final 
periodic attractor following the given initial state. The 
timescale of these patterns becomes large for networks 
with high K, making it more difficult to keep the exact 
dynamic pattern under the mutation of a weight. In pop- 
ular terms, an increased complexity of the network makes 
further evolution difficult. One may speculate that this 
is the reason for real genetic networks to keep their con- 
nectivity low: It will be easier to evolve by increasing the 
number of genes N at a fairly low connectivity level (the 
present model, however, does not consider variable N). 

In Fig. la we further see that marked punctuations 
occur, where long periods of nearly fixed average connec- 
tivity sometimes are interrupted by a sudden change in 
connectivity. This interplay between long waiting times 
and short times for actual changes is in fact observed 
in the fossil record. The phenomenon has been coined 
"punctuated equilibrium" by Gould and Eldredge |J . As 
also seen from Fig. lb, the periods of stasis show a sim- 
ilar structure on shorter timescales as they do on longer 
timescales. This is explored further in Fig. 3a where we 
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FIG. 3. Stasis time distribution in the neutral evolution 
of networks (a) and the decomposition into stasis time distri- 
butions for different intervals of the average connectivity (b). 
The simulation is the same as in Figure 1. While the aver- 
age over all networks approximately follows a power law, the 
distributions for networks with restricted K do not. One ob- 
serves that large connectivity typically implies a lower degree 
of evolvability. 

Approximately the distribution of stasis times is oc 
1/t 2 . Periods of stasis at high values of K can become 
long, which in practice calls for very long equilibration 
times. 

In Fig. 3b we decompose the stasis time distribution 
into times obtained for different values of the average 
connectivity. Again we observe that higher K typically 
shows longer stasis times. Remarkably, when looking at 
the statistics of a small interval of low K values, we ob- 
serve exponentially distributed stasis times. The power 
law behavior then comes about by averaging over the 
range of all K values. 

In order to test for the robustness of our model we 
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tried other mutation rules (again without any evolution- 
ary pressure, i.e., symmetric in adding and removing 
weights) . In one variant a daughter network was created 
by a) adding or b) removing a weight in the coupling ma- 
trix at random, with p = 1/2 each, thus allowing for K- 
changing mutations only. In another variant a daughter 
network was created by independently adding a random 
weight with p = 1/2 and removing a random weight with 
p = 1/2. We also tested a scenario where we demanded 
complete match between attractors of mother and daugh- 
ter for two different initial configurations. Also, we con- 
sidered this case with demanding only partial overlap be- 
tween mother and daughter, i.e., match in at least one 
of the two tested environments, only. Finally, we tested 
networks with weights Wij E { — 1,0,1} such that the sig- 
nal transmitted from an inactive node can differ from 
the value of a disconnected input node. In all cases our 
results were robust. 

Let us briefly discuss the meaning of the stasis times 
and punctuations observed here. According to the defini- 
tion of our model, we quantify the waiting time in terms 
of the number of times mutant networks are exposed to 
new environments before a neutral mutation occurs that 
fulfills continuity. Thus they are not to be confused with 
the "neutral evolution" introduced by Kimura |$) which 
leads to waiting times consisting of a number of neutral 
mutations. The genetic networks are formally defining a 
"species" and the length of the waiting times indicates 
the "genetic flexibility" of a species. 

Associating the interconnectedness of the networks 
with the genetic flexibility of real organisms one may 
attempt to understand a puzzling decomposition of life- 
times of species in the fossil record. First it was noted by 
Van Valen |l6| that each group of closely related species 
has exponentially distributed lifetimes. Second, an anal- 
ysis of the overall distribution of genera lifetimes, tab- 
ulated by Raup and Sepkos ki Il7[ , showed that this is 
rather distributed as oc 1/t 2 |l9| for genera lifetimes ex- 
ceeding 10 million years. It is tempting to speculate 
that groups of closely related species are associated to 
the same genetic flexibility, and thus evolve, and eventu- 
ally get extinct, with a frequency given by this genetic 
flexibility. This would explain the exponential distribu- 
tion of Van Valens. Averaging over all genetic flexibilities 
is then an average over different characteristic lifetimes, 
and our simplified evolution scenario demonstrates how 
such an averaging can give an overall 1 /t 2 distribution. 

The obtained 1/t 2 scaling may be an inherent part of 
our neutral evolution scenario jlq] . In comparison, for 
evolution on fitness landscapes one typically obtains a 
distribution oc 1/t corresponding to a sampling of wait- 
ing times for passing over barriers ]j"9pc| ] , although a 
1 /t 2 distribution can be obtained by supplementing a hill 
climbing concept with the assumption that extinction of a 
given species is determined by evolutionary rates of older 
species pl[. Abandoning fitness landscapes, the ecologi- 



cal network model of ref. [|22j|23|] instead determines the 
fate of a species through an ecological connectivity ma- 
trix. When this evolving network of species is assigned 
a connectivity that is comparable to (eco)system size, it 
shows a 1/t 2 distribution of genera lifetimes. In con- 
trast to these macro-evolutionary models, our study of 
Boolean networks only considers one species, with a com- 
parison to species extinction data that is based on an 
extrapolation of the obtained evolutionary clock. 

In conclusion we have studied evolution of Boolean net- 
works in absence of any competition. This simplification 
allowed us to discuss how the requirement of evolving ro- 
bust networks in itself may lead to an evolution which 
exhibits punctuated equilibrium. 
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